با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید. 

library(plotly)
library(readr)
library(gganimate)
library(ggmap)
library(dplyr)
library(ggplot2)
library(highcharter)
library(ISOcodes)
library(countrycode)
library(stringr)
library(lubridate)
equake = read_csv("data/worldwide.csv")
Country = c()
for(i in 1:length(equake$place)){
  Country =c(Country, tail(str_split(equake$place[i], ", ")[[1]], 1))
}
equake$Country = Country
equake$Year = format(equake$time, "%Y")
equake$Month = format(equake$time, "%m")
equake$Day = format(equake$time, "%d")

iequake = read_rds("data/iran_earthquake.rds")
disaster = read_delim("data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
sequake = read_rds("data/historical_web_data_26112015.rds")

۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید. 

plot_ly(sequake, x=~Longitude, y=~Latitude, z=~Depth, size=~Magnitude)
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
## Warning: package 'bindrcpp' was built under R version 3.4.4

۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.) 

worldMap <- map_data("world")
animateData = disaster %>% filter(FLAG_TSUNAMI == "Tsu" & !is.na(LATITUDE) & !is.na(LONGITUDE) & !is.na(LOCATION_NAME)& !is.na(EQ_PRIMARY)) %>% arrange(EQ_PRIMARY)
p = ggplot() +
    geom_polygon(data = worldMap, aes(long, lat, group=group), fill = "white", color = "lightblue") +
  geom_point(data = animateData, aes(LONGITUDE, LATITUDE,
                                                  frame = paste0("Magnitude under ", floor(EQ_PRIMARY)),
                                                  cumulative = T, size = EQ_PRIMARY, color = EQ_PRIMARY)) + scale_color_continuous(low = "yellow", high = "red", guide = F) + scale_size(guide = F)
gganimate(p, filename = "tsunami.gif")

۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید). 

iranMap = read_rds("data/Tehrn_map_6.rds")
ggmap(iranMap) + stat_density_2d(geom = "polygon", data = iequake, aes(x =Long, y = Lat, alpha = ..level.., fill = ..level..)) + scale_fill_gradient(low = "blue", high = "red", guide=FALSE) +  scale_alpha(range = c(0.5, 1), guide = FALSE)

۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)

برای این سوال بنده روش های مختلفی را بررسی کردم و در نهایت به این نتیجه رسیدم اینکار می تواند کار منطقی ای باشد.
به اینصورت که توزیع تجربی فواصل میان زلزله های مختلف بالای ۷ ریشتر که در ایران به وقوع پیوسته را در نظر بگیرم و بر اساس آن احتمال بودن زلزله بعدی تا ۵ سال آینده را از تابع احتمال تجمعی بیابم. 

iequake$Year = format(iequake$OriginTime, "%Y")
yearDiffs = disaster %>% filter(COUNTRY == "IRAN" & EQ_PRIMARY > 7) %>% mutate(previousYear = lag(YEAR)) %>% filter(!is.na(previousYear)) %>% mutate(diff = YEAR - previousYear) %>% .$diff
empiricalCDF = ecdf(yearDiffs)

با توجه به اینکه آخرین زلزله در سال ۲۰۱۷ رخ داد، برای محاسبه احتمال رخ دادن زلزله بزرگ تا ۵ سال آینده باید احتمال کمتر بودن فاصله زلزله بعد از ۶ و بیشتر بودن آن از ۱ را بیابیم. پس احتمال آمدن زلزله در پنج سال آینده به شرط نیامدن زلزله در یک سال اخیر را می یابیم: 

(empiricalCDF(6) - empiricalCDF(1))/(1-empiricalCDF(1))
## [1] 0.173913

۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!) 

deaths = disaster %>% filter(!is.na(TOTAL_DEATHS) & !is.na(COUNTRY)) %>% group_by(COUNTRY) %>% summarise(TotalDeaths = sum(TOTAL_DEATHS), AverageDeaths = mean(TOTAL_DEATHS)) %>% select(country = COUNTRY, TotalDeaths, AverageDeaths)

deaths$Code = countrycode(deaths$country, "country.name", "iso3c")

hcmap(data = deaths, joinBy = c("iso-a3", "Code"), value = "TotalDeaths", name = "Total Deaths")
hcmap(data = deaths, joinBy = c("iso-a3", "Code"), value = "AverageDeaths", name = "Average Deaths")

۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید. 

trainData = disaster %>% filter(!is.na(TOTAL_DEATHS) & !is.na(LATITUDE) & !is.na(LONGITUDE) & !is.na(FOCAL_DEPTH) & !is.na(EQ_PRIMARY))
fit = lm(data = trainData, TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY)
summary(fit)
## 
## Call:
## lm(formula = TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + 
##     EQ_PRIMARY, data = trainData)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -14162  -4354  -2228     19 311179 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.837e+04  3.599e+03  -5.105 3.93e-07 ***
## LONGITUDE   -2.355e-01  6.757e+00  -0.035  0.97221    
## LATITUDE     6.956e+01  2.593e+01   2.682  0.00743 ** 
## FOCAL_DEPTH -2.379e+01  1.384e+01  -1.719  0.08596 .  
## EQ_PRIMARY   3.172e+03  5.370e+02   5.907 4.72e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17780 on 1047 degrees of freedom
## Multiple R-squared:  0.03526,    Adjusted R-squared:  0.03157 
## F-statistic: 9.567 on 4 and 1047 DF,  p-value: 1.342e-07

۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟

برای اینکار ابتدا تعریفی برای پیش لرزه ارائه می دهیم و سپس سعی می کنیم به کمک پیش لرزه ها و مدل رگرسیون خطی، پیش بینی ای برای شدت و زمان زلزله اصلی ارائه کنیم و دقت آن را بررسی کنیم. پیش لرزه را اینگونه تعریف می کنیم که برای هر کشور در هر روز، بزرگترین زلزله را به عنوان زلزله اصلی و سایر زلزله ها را به عنوان پیش لرزه و پس لرزه در نظر می گیریم. سپس زلزله های غیر اصلی که قبل از زلزله اصلی آمدند را به عنوان پیش لرزه برای فیت کردن مدل نگه داری می کنیم. 

foreShockData = equake %>% group_by(Country, Year, Month, Day) %>% mutate(isMain = (mag == max(mag)), mainMag = max(mag)) %>% arrange(Country, Year, Month, Day, desc(isMain)) %>% mutate(mainTime = first(time)) %>% mutate(timeDiff = floor(as.integer(mainTime - time)/3600)) %>% filter(mainTime>time)

برای هر پیش لرزه شدت زلزله اصلی و فاصله زمانی ساعتی تا وقوع زلزله اصلی را نگه داری می کنیم و سعی می کنیم این مقادیر را با کمک شدت و عمق زلزله حدس بزنیم. 

#Shuffling data
foreShockData = foreShockData[sample(nrow(foreShockData)),]
train = foreShockData[1:as.integer(0.9 * nrow(foreShockData)),]
test = foreShockData[-(1:as.integer(0.9 * nrow(foreShockData))),]
fit = lm(data = train, mainMag ~ depth + mag)
summary(fit)
## 
## Call:
## lm(formula = mainMag ~ depth + mag, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2212 -0.4577 -0.1280  0.3619  4.0278 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.9923229  0.0264461  75.335  < 2e-16 ***
## depth       -0.0003755  0.0000588  -6.386 1.76e-10 ***
## mag          0.6983743  0.0074985  93.136  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6726 on 13308 degrees of freedom
## Multiple R-squared:  0.4027, Adjusted R-squared:  0.4026 
## F-statistic:  4485 on 2 and 13308 DF,  p-value: < 2.2e-16
summary(stats::predict(fit, test)-test$mainMag)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -3.384082 -0.339906  0.127882 -0.004059  0.461447  1.195886
timeFit = lm(data = train, timeDiff ~ depth + mag)
summary(timeFit)
## 
## Call:
## lm(formula = timeDiff ~ depth + mag, data = train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.460 -5.040 -1.187  4.023 16.237 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.8604321  0.2273290  34.577  < 2e-16 ***
## depth        0.0007553  0.0005054   1.494 0.135075    
## mag         -0.2255332  0.0644564  -3.499 0.000469 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.782 on 13308 degrees of freedom
## Multiple R-squared:  0.0009507,  Adjusted R-squared:  0.0008005 
## F-statistic: 6.332 on 2 and 13308 DF,  p-value: 0.001784
summary(stats::predict(timeFit, test)-as.integer(test$timeDiff))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -16.0577  -3.7726   1.2509   0.1406   4.9405   7.4125

مدل را بر روی ۹۰ درصد داده فیت می کنیم و به کمک ۱۰ درصد دیگر خطا ها را بررسی می کنیم. مدل مربوط به شدت زلزله اصلی R-squared حدود ۰.۴۰ و مدل مربوط به فاصله زمانی R-squared بسیار پایینی دارد. این نشان دهنده این است که با اینکه هر دو مدل مدل های نادقیقی برای پیش بینی هستند اما پیش بینی کننده شدت خطای بسیار کمتری نسبت به پیش بینی کننده زمان دارد.
پس از انجام پیش بینی بر روی داده تست در مدل شدت می بینیم به طور میانگین پیش بینی ها حدود ۰.۰۰۷ ریشتر بیشتر از مقدار واقعی هستند بیشتر خطاها بین مثبت منفی ۰.۵ ریشتر هستند. این خطا در مدل مربوط به زمان به گونه ای است که به طور میانگین خطای پیش بینی حدود ۰.۰۰۲ ساعت است اما بیشتر خطاها بین مثبت منفی ۴.۵ ساعت هستند که خطای قابل توجهی می باشد!

۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.) 

data = equake %>% filter(!is.na(depth) & !is.na(mag))
cor.test(data$depth, data$mag)
## 
##  Pearson's product-moment correlation
## 
## data:  data$depth and data$mag
## t = 50.03, df = 60629, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1914584 0.2067469
## sample estimates:
##       cor 
## 0.1991147
hchart(equake %>% head(1000), type = "point", hcaes(depth, mag)) %>% 
  hc_add_theme(hc_theme_darkunica())

با توجه پایین بودن مقدار کوریلیشن به دست آمده، نتیجه می گیریم همبستگی بین این دو متغیر پایین است و شدت زلزله به عمق آن بستگی ندارد. نمودار نقطه ای نشان دهنده شدت زلزله بر حسب عمق زلزله نیز گویای این مسئله می باشد.

۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد. 

equake %>% group_by(Country, Year) %>% summarise(AverageMagnitude = mean(mag, na.rm = T))
## # A tibble: 925 x 3
## # Groups:   Country [?]
##    Country     Year  AverageMagnitude
##    <chr>       <chr>            <dbl>
##  1 Afghanistan 2015              4.33
##  2 Afghanistan 2016              4.35
##  3 Afghanistan 2017              4.30
##  4 Afghanistan 2018              4.35
##  5 Alabama     2016              2.75
##  6 Alabama     2017              2.66
##  7 Alaska      2015              2.99
##  8 Alaska      2016              3.08
##  9 Alaska      2017              3.11
## 10 Alaska      2018              3.18
## # ... with 915 more rows

از آن جایی که داده های مربوط به سال های اخیر در اختیار ماست و پروژه هارپ پروژه ای قدیمی است نمی توانیم به طور دقیق آن را بررسی کنیم. اما با توجه به همین داده ها و با این فرض که آمریکا قصد ایجاد تاثیر بر روی برخی کشورهای دیگر دارد، با کمک گرفتن از ایران به عنوان یکی از کشورهای مشکوک، نحوه تغییر تعداد زلزله ها و شدت آن ها را برای سال های مختلف در نظر می گیریم و با آمریکا مقایسه می کنیم. 

iran_equake = equake %>% filter(Country == "Iran")
us_states = "Alabama\nAlaska\nArizona\nArkansas\nCalifornia\nColorado\nConnecticut\nDelaware\nFlorida\nGeorgia\nHawaii\nIdaho\nIllinois\nIndiana\nIowa\nKansas\nKentucky\nLouisiana\nMaine\nMaryland\nMassachusetts\nMichigan\nMinnesota\nMississippi\nMissouri\nMontana\nNebraska\nNevada\nNew Hampshire\nNew Jersey\nNew Mexico\nNew York\nNorth Carolina\nNorth Dakota\nOhio\nOklahoma\nOregon\nPennsylvania\nRhode Island\nSouth Carolina\nSouth Dakota\nTennessee\nTexas\nUtah\nVermont\nVirginia\nWashington\nWest Virginia\nWisconsin\nWyoming\nDistrict of Columbia\nPuerto Rico\nGuam\nAmerican Samoa\nU.S. Virgin Islands\nNorthern Mariana Islands"
us_equake = equake %>% filter(str_detect(us_states,Country))
iran_us = rbind(iran_equake, us_equake %>% mutate(Country = "United States"))

#print chart
hcboxplot(x = iran_us$mag, var = iran_us$Year, var2 = iran_us$Country) %>% hc_chart(type = "column")
iran_us %>% group_by(Country, Year) %>% summarise(count = n()) %>% hchart("line", hcaes(x = Year, y = count, group = Country))

با توجه به نمودار جعبه ای مربوط به شدت زلزله ها و نمودار مربوط به تعداد زلزله های سالانه افزایش غیرطبیعی ای در ایران مشاهده نمی شود و به نظر نمی توان این فرضیه را تایید کرد.

۱۰. سه حقیقت جالب در مورد زلزله بیابید. 

equake$hour = format(equake$time, "%H")
equake %>% group_by(hour) %>% summarise(averageMag = mean(mag), count = n()) %>% hchart("point", hcaes(x = "hour", y = "count", size = "averageMag"), color = "black")
hcboxplot(x = equake$mag, var = equake$hour, color = "black")

در نمودار بالا زلزله های جهان به تفکیک ساعت بررسی شده اند تا تعداد و شدت آن ها مقایسه شوند.
با توجه به نمودارهای به دست آمده می توان دید بین ساعت ۱۲ تا ۱ بامداد از خطرناک ترین زمان هاست که در آن علاوه بر بالا بودن تعداد زلزله ها، زلزله ها از شدت بالایی نیز برخوردار بوده اند.
در مقابل با توجه به نمودار تعداد زلرله بر حسب ساعت و توزیع شدت زلزله ها بر حسب ساعت به نظر می رسد در ساعات ۶ تا ۷ صبح نیز هم خطر وقوع زلزله و هم خطر بالا بودن شدت آن کمتر از ساعات دیگر است.
در ساعات دیگر نیز می توان دید تقریبا توزیع شدت زلزله ها ثابت است و با رفتن از ساعات اولیه روز به ساعات انتهایی تعداد زلزله ها افزایش می یابد. 

iequake$hour = format(iequake$OriginTime, "%H")
iequake %>% group_by(hour) %>% summarise(averageMag = mean(Mag), count = n()) %>% hchart("point", hcaes(x = "hour", y = "count", size = "averageMag"), color = "green")
hcboxplot(x = iequake$Mag, var = iequake$hour, color = "green")

با انجام همین بررسی در زلزله های ثبت شده در ایران در سال ۲۰۱۵، نتایج به دست آمده را بررسی می کنیم.
با توجه به نمودار تعداد بر حسب ساعت می بینیم که در ساعات ۶ تا ۱۵ شاهد یک قله در تعداد هستیم که نشان می دهد بیشترین تعداد زلزله ها در ساعات ۹ تا ۱۰ صبح به وقوع می پیوندند.
میانگین شدت در ساعات مختلف نیز تقریبا به هم نزدیکند ولی این میزان در ساعات ۹ تا ۱۰ کمتر و پراکندگی آن نیز کمتر است. 

iequake %>% group_by(hour) %>% summarise(averageMag = mean(Mag), count = n()) %>% hchart("point", hcaes(x = "count", y = "averageMag",  group = "hour"), color = "green") %>% hc_legend(enabled = F)
iequake %>% group_by(hour) %>% summarise(averageMag = mean(Mag), count = n()) %>% cor.test(formula = ~ averageMag + count, data = .,  method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  averageMag and count
## S = 4416, p-value = 2.489e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##   rho 
## -0.92

با توجه به نتایج به دست آمده از زلزله های ایران حقیقت جالبی به نظر رسید که آن را بررسی می کنیم.
به نظر می رسد در ایران با افزایش تعداد زلزله ها در ساعات مختلف، شدت آن ها کاهش می یابد.
این فرضیه با توجه به نمودار میانگین شدت بر حسب تعداد و تست کوریلیشن رنک شده بر روی این دو متغیر تایید می شوند!